clear; clc;

m = 100;
n = 50;
maxIter = 10000;

A = randn(m, n);
k = 2;

[ U, S, V ] = svd(A, 'econ');
S = diag(S);
sk1 = S(k + 1);
% S(1:k+1) = sk1;
A = U*diag(S)*V';
[ U, S, V ] = svd(A, 'econ');
U = U(:,1:k);
S = diag(S);

R = rand(n, k);
Y = A*R;
[Q, ~] = qr(Y, 0);
I = eye(size(Q, 1));
fberr = zeros(maxIter, 1);
sperr = zeros(maxIter, 1);
for i = 1:maxIter
    fberr(i) = norm(Q*Q'*A - U*U'*A, 'fro');
    sperr(i) = norm(Q*Q'*A, 'fro') - norm(U*U'*A, 'fro');
    
    Y = A*(A'*Q);
    [Q, ~] = qr(Y, 0);
    fprintf('err %d \n', fberr(i));
    if(i > 1 && abs(fberr(i) - fberr(i - 1)) < 1e-10)
        break;
    end
end

clear Xapp Xopt;

sperr = sperr(1:i);
fberr = fberr(1:i);

close all;
figure;
semilogy(abs(fberr));
figure;
semilogy(abs(sperr));

[U, ~, ~] = svd(A, 'econ');
U = U(:, 1:k+1);